
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!the funtion of dot product of two vector
function prod(s1,s2) result(s3)
implicit none
double precision,dimension(3)::s1,s2
double precision::s3
integer::i
s3=0
!write(*,*)s1,s2
do i=1,3
s3=s3+s1(i)*s2(i)
end do
end function
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!calcaulate initial magnitization
subroutine initial_mag
use para
implicit none

integer::i,j
totmx=0.0
totmy=0.0
totmz=0.0
do i=1,nlatt
  do j=1,natom
  totmx=totmx+spin(1,i,j)
  totmy=totmy+spin(2,i,j)
  totmz=totmz+spin(3,i,j)
  end do
end do
  totm2=totmx*totmx+totmy*totmy+totmz*totmz
  totm=sqrt(totm2)
end subroutine


!calculate the initial energy
subroutine initial_energy
use para
implicit none
double precision,external::prod
integer::i,j,k,l1,m,n
real(kind=8)::sumjexc

tote = 0.0
dote = 0.0
sumjexc=0.0
do i = 1, nlatt, 1
  do j=1, natom,1
    do k=1, natom,1
      do l1=-mmldx, mmldx,1
        do m=-mmldy,mmldy,1
          do n=-mmldz,mmldz,1
            if(n==0 .and. &
              &m==0 .and. &
              &l1==0 .and. &
              &j==k) then
             if(my_rank==0 .and. i==1) write(*, "(5I3, F10.2)") &
             &j,k,n,m,l1,jexc(j,k,n,m,l1)
            else
            dote(j)=dote(j)+&
            &jexc(j,k,n,m,l1)*&
            &(prod(spin(1,i,j), &
            &spin(:,&
            &mld(n,m,l1,i),k)))
            if(my_rank==0) sumjexc=sumjexc+jexc(j,k,n,m,l1)
            end if
          end do ! n
        end do ! m
      end do ! l1
    end do ! k
  end do !j
end do ! i

do i=1, natom, 1
  tote = tote + dote(i)
  if(my_rank==0) write(*,"(A7, I3, F12.3)") "DOTE",i,dote(i)
end do
tote=tote/2.0
sumjexc = sumjexc / 2.0

if(my_rank==0) then
 write(*,"(A7,F12.3,A10,F12.3)") "TOTE:",tote,"SUMJEXC:",sumjexc
 write(*, "(A9,2X,3F6.2)") "SPIN1:", spin(:, 1, 1)
 write(*, "(A9,2X,3F6.2)") "SPIN2:", spin(:, 1, 2)
 write(*, "(A9,2X,3F6.2)") "SPIN3:", spin(:, 1, 3)
 write(*, "(A9,2X,3F6.2)") "SPIN4:", spin(:, 1, 4)
 write(*, "(A9,2X,3F6.2)") "SPIN5:", spin(:, 1, 5)
 write(*, "(A9,2X,3F6.2)") "SPIN6:", spin(:, 1, 6)
end if
!############################# add the anisotropy energy
if(anis_flag==1)then

do i=1,nlatt
   do j=1,natom
   tote=tote-(prod(D4,spin(1,i,j)))**2
   end do
end do

end if

if(hflag==1)then

do i=1,nlatt
  do j=1,natom

  tote=tote+prod(hfield,spin(1,i,j))*magmom

  end do
end do

end if

end subroutine


subroutine initial_count
use para
implicit none
 uacceptb=0
 ucountb=0
 uacceptl=0
 ucountl=0
 scount=0
 sacpt=0
 
end subroutine
